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We analytically study a system of spinless fermions driven at the boundary with an oscillating 
chemical potential. Various transport regimes can be observed: at zero driving frequency the particle 
current through the system is independent of the system's length; at the phase-transition frequency, 
being equal to the bandwidth, the current decays as ~ n~ a with the chain length n, a being either 
2 or 3; below the transition the scaling of the current is ~ n" 1 ' 2 , indicating anomalous transport, 
while it is exponentially small ~ exp (— £•) above the transition. Therefore, by a simple change of 
frequency of the a.c. driving one can vary transport from ballistic, anomalous, to insulating. 
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PACS numbers: 05.60.Gg, 03.65.Yz, 75.10.Pq 



I. INTRODUCTION 

A system of free fermionic particles is of great impor- 
tance in many areas of physics. It can describe a gas of 
electrons, a solid state system in a tight-binding approx- 
imation, or the so-called XX spin chain, to name just 
a few. In a noninteracting case without fields physics is 
rather simple. Free modes, either in real or in momentum 
space, propagate without mutual interaction or dissipa- 
tion with constant speed, resulting in ballistic transport 
properties. Things become more interesting in the pres- 
ence of external fields. Limiting the discussion to one- 
dimensional (ID) systems, any on-site disorder results 
in a localization of all states - the Anderson localiza- 
tion [I]. A constant static electric field causes Bloch os- 
cillations [2] in a periodic lattice. An oscillating electric 
field on the other hand results in a renormalized hopping 
strength, leading to localization at specific values of the 
field or the frequency [3]. In the present work we shall 
study a new setting, namely that of a harmonically (a.c.) 
boundary-driven system of free fermions, and show that, 
depending on the driving frequency, very different trans- 
port regimes can be observed. The model can also be 
seen as a minimal model for the boundary driven non- 
equilibrium quantum phase transition. We also discuss 
its experimental implementations, showing that it could 
be realized in an optical setting with cold atoms or ions 
or, in its fermionic version, with mesoscopic systems. 

We shall study a ID chain of spin-1/2 particles inter- 
acting between nearest-neighbor sites with an exchange 
interaction of the XX type. The Hamiltonian is 



coupled to magnetization reservoirs. In the language of 
fermions the reservoirs impose an external chemical po- 
tential determining the number of fermions in the system. 
The evolution of the system's density matrix is described 
by the Lindblad master equation pQ, 

dp/dt = i[p,H]+C dis (p), (2) 

where a dissipator £ dls (p) = Y,k {[ L kP,L k ] + [L k ,pL k ]) 

takes into account the influence of two reservoirs. The 
left reservoir is described by a pair of time-dependent 
Lindblad operators Li. 2 (t) = y/e{l ± pjt))af, and sim- 
ilarly at the right end, L 3 ^(t) — y/e(l T /•*(*)) &„. Note 
that drivings at the two ends are opposite. For instance, 
when there is a maximal coefficient of a^ at the left 
end, we have a minimal coefficient of <7+ at the right 
end. Two important bath parameters arc the coupling 
strength e between the chain and the reservoirs and the 
value of the "chemical potential" that is oscillating in 
time, p = pocosojt. Nonzero p means that probabilities 
for an injection and absorption of a spin/fermion are dif- 
ferent. For d.c. driving, u> = 0, the transport of magne- 
tization is trivially ballistic, i.e., the value of the current 
through the system is independent of the length n, and 
the nonequilibrium stationary solution of the Lindblad 
equation is known [5] and can be compactly written in 
a matrix product operator form [5J. We shall show by 
a fully analytic calculation that the situation is rather 
different and interesting for uj ^ 0. 



II. THE SOLUTION 
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Equivalently, using Jordan- Wigner transformation the 
XX model can be mapped to a system of non- 
interacting spinless fermions with the Hamiltonian H = 
'•' ' where c k = ^kHj<k a p and 
(<x£ ± i(j y k )/2) satisfy the stan- 
The first and the last spin are 
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dard fermionic algebra. 



We want to find an asymptotic solution p(t — > oo) 
of the Lindblad equation to which the system converges 
after very long time. Because the system is harmonically 
driven this solution is not time- independent, as is the 
case for u) — 0, but instead oscillates with the very same 
forcing frequency u>. In fact, the asymptotic stationary 
solution of Eq. pi can be written in the form 



p(t) = 2-"l 



\{P^ 
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where p is time-independent. The nonequilibrium sta- 
tionary state p is in our case unique. There are sev- 
eral ways of computing p. One is to realize that the 
XX model, even when time-dependent, belongs to a 
class of systems [?H5] in which exponentially many equa- 
tions for all r-point functions decouple into a hierar- 
chy in which one has separate equations for each or- 
der of correlations. This greatly simplifies the Lindblad 
equation, enabling one to find the exact non-equilibrium 
steady state [8] and even study time-evolution towards 
a steady state [5]. We shall be especially interested in 
2-point correlations as these include local magnetization 
a k (=fcrmion density) as well as spin current (^particle 
current) j k = 2(a k a k+1 — & k cr k+1 )- In fact, the only 

nonzero 2-point terms in p are B k = 2{c k c k + r — c k c k+r ) 
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for even r and B k — 2i(c k c k + r + CfcCJL r ) for odd r 
[IB] , For instance, the first two operators in the series 
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are B k ' = — 2a k and B k ' = j k /2. In a chain of length 
n there are in total n(n — l)/2 such operators and one 
can write a self-contained set of as many equations for 
their unknown coefficients in p. Due to our parametriza- 
tion of driving they are all exactly proportional to po; 
from now on we therefore set /zo = 1 in all our results. 
Using standard procedures such a set of linear equations 
can be solved for chain lengths of n ~ 10 3 . Once these 
2-point coefficients are known, they can be used as in- 
homogeneous source terms in the set of equations for all 
4-point terms (which are all proportional to pfy, and so 
on (see [5] for an example of such a calculation). 

In our case though one can do even better. Follow- 
ing [T0UT2"] the evolution of 2-point correlations (i.e., 
the above mentioned n(n ~ l)/2 linear equations) can 
be compactly written in terms of a linear matrix equa- 
tion. Defining the time dependent covariances (wjWk) = 
tr p(t)WjW k = 6 jt k-iZj,k, where w m = c m + cl nl w m+n = 
i( c m — c ln)i Tn = 1, . . .n, and applying the Lindblad mas- 
ter equation, we obtain a differential equation for the 
2n x 2n covariance matrix Z 



dZ/dt 



-X x Z-ZX + Ycoswt, 



(4) 



where X = 2i u y ® J - 2el 2 <8> R and Y = -4i a y ® P are 
real matrices, with J,P,R being n x n matrices having 
nonzero elements Jk,k+i = Jk+i,k = 1, fc = 1, • • • ,n — 1, 
and Ria = i?„,„ = P^i = -P„,„ = 1. By consid- 
ering the symmetry of the differential equation Q we 
write the covariance matrix in the simple form Z(£) = 
Re (e i "'(l 2 <8> Z - i a y ® Z 2 )). Inserting the last defini- 
tion into Q we get two coupled equations for the n x n 
matrices Zo and Z 2 . We immediately see that (Zo)j.fc = 
for j + k even, and (Z 2 )j ! fc — for j + k odd, hence we 
can replace the two equations by one. Defining the or- 
thogonal transformation Oj : k = (—^) k+1 Sj } k, the matri- 
ces C* = 0(Zo ± iZ 2 ), and inverting the last relation, 
we end up with two uncoupled equations 

2{J±rcR,C =F }=FwC =F = -4eOP, (5) 

where {, } is an anticommutator. Note that for odd n, 



OP = P, while for even n, OP = R. The last equation 
([5]) can be solved perturbatively in the coupling e. The 
solution to first order in the coupling, obtained by the 
Fourier method (see the Appendix [b]) , is 



c 7. k = 
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p,m— 1 
p+m— n(mod 2) 



sin a p sin a rn sin a,j p sin a^ 
(n+l) 2 (X p + X m ) 



(6) 



( | - 4 cos a m ) - ^j sin a m where a k = ^ . 
It is instructive to emphasize the connection Cj h — 
-i(-iy {wjW k ) = -|(-1) J '{B$ k ~ jl) ) iij + k is odd and 

Cl k = (-iy(w J+n w k ) = -i(-l)i<Bf-^)) iij + k is 
even. All expectation values in this work are computed 
with respect to the a.c. part of the density matrix p, 



(A) = tT(pA). Note that C+ k = {-l) J+k+1 Cj k and 
that n x n matrix C _ contains all non-vanishing matrix 
elements of Z. 

From the perturbative solution ffib we can see [14] that 
for large n and uj > lo c = 8 the imaginary part of the Lya- 
punov equation ([5]) can be neglected. In the thermody- 
namic limit one can in fact express the solution in terms 
of a generalized hypergeometric function (Appendix O) . 
Here we prefer a different route by first taking a con- 
tinuum limit and writing a partial differential equation. 
This results in a Helmholtz equation for the correlation 
function C(x = -^rr,y = :ttt) = C, : 



V 2 C(x, y) + (« c 



-^-1 - V 

n+1' j,k> 

oj)n 2 C(x,y) 



P(x,y), (7) 



whexeP(x,y)=e(6(x-±,y-±) + (-l) n 8(x-l + ±,y- 

1 + -)) is the source term resulting from the driving. 
Since the dissipation effects in the matrix X can be ne- 
glected, the boundary conditions arc C(0, y) = C(x, 0) = 
C(l,y) = C(x, 1) = 0. Using the method of mirror im- 
ages the problem can be unfolded into an infinite plane, 
where the source term becomes a grid of quadrupolcs 



The solution is 



C(x,y) =e 



CSO 



(-\Y n 
/ ,^—G{x-2j-v,y-2k-v), (8) 
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where G(x, y) is the Green's function of Eq. ([7]) for a 
single quadrupole source q{x,y) = 25' (x)S' (y) at the ori- 
gin, where S'(x) denotes the derivative of Dirac's function 
in the sense of distribution. For the critical frequency 
lu = uj c the Green's function is particularly simple (and 
ri- independent) 



G c (x,y)=4xy/[Tr(x 2 + y 2 ) 2 



(9) 



For u> < w c , we can use heuristic arguments, combining 
integral representation of the sum (pSJ) and stationary- 
phase integration leading to a universal asymptotic 
n— scaling of the covariances 



\C~ k \ = Oin- 1 ' 2 ), if i,fecxn. 



(10) 



III. RESULTS 

The main quantity we shall be interested in is the cur- 
rent in the middle of the chain, Ji ( n +i)/2J > an d its scaling 
with n. The scaling will tell us how much the disturbance 
from the driving at the chain ends influences the bulk. In 
Fig. [I] we show the current's dependence on the driving 
frequency uj, which turns out to be the most important 
parameter in the system, that can qualitatively change 
the transport. For small couplings e, shown in the top 
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FIG. 1. Dependence of the current in the middle of the chain 
on the driving frequency uj. Top: For small coupling s = 
0.0I,n = 257, many resonances are visible. Bottom: For 
strong coupling e = I we have a smooth dependence. For 
uj > uj c (lower inset) the current decays exponentially with n. 

frame of Fig. [I] , one can see that overall the current de- 
cays with uj, however for uj < uj c , there are many very 
narrow resonances. Looking at the theoretical formula 
ffity, valid for small s, it is easy to understand where these 
resonances come from. They arise from points where 
the denominator A p + A m is very small, which happens 
when uj — uj p - m = e p + e m , with e m = 4cosa m being 
the energies of free fermionic modes. The width of res- 
onances scales as Aw ~ e/n, and within the resonance 
region \uj — uj p ^ m \ < Aw the scaling of covariances is 



Cj k \ ---- O^n- 1 ) rather than ©(e 1 ^ 1 / 2 ) |f0|). For 



larger couplings (bottom frame in the Fig. [I]) the reso- 
nances merge into a smooth curve. We can also see that 
for uj > uj c the current decays very rapidly with uj and n, 
so there the system becomes an insulator for n — > oo. At 
uj = uj c we have a non-equilibrium phase transition. The 
value of uj c = 8 is also given by the bandwidth of the 
XX chain or, in physical picture, an insulating regime 
appears because fermions do not have enough time to 



reach the neighboring site before the sign of the driv- 
ing reverses. Next, we show in Fig. [2] the spatial de- 
pendence of the magnetization (cr|) and the spin current 
(jk) along the chain. Note that p(t) is time-dependent, 
even after a long time, so the continuity equation reads 
kj(er|) = (jk-i) — (jk)- In addition to the amplitudes, 
shown in Fig. [21 there is also a nontrivial dependence of 
phases on the spatial position k of each of these quanti- 
ties. We should note that none of the qualitative features 
depends on the value of e (for more see the Appendix [A|. 
At the transition point uj — uj c the correlations can be for 
large n approximately calculated by summing over only 
two nearest quadrupoles in dsl p)| , resulting in 



C c (x,y) 



(G c (x,y) + (-l) n G c (x-l,y-l)). (11) 



The difference between this theoretical dependence and 
the exact values is shown in Fig. [2] and is small. For 
uj < uj c the main wavelength visible in magnetization and 
current profiles is determined by the closest resonance. 
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FIG. 2. (Color online) Magnetization profiles (left column) 
and spin current profiles (right column) for n = 257, e = 0.1, 
and three different uj. At the critical uj c = 8, the profiles agree 
with the theoretical estimate (III (relative error is always 
less than 10%). Dashed line for uj — 8.1 is exp (— fc/£) with 

To assess the nature of spin transport at various val- 
ues of uj we have studied dependence of the spin current 
in the middle of the chain on the system size n. Above 
the transition point uj > uj c the dependence is exponen- 
tial, |i(„ + i)/ 2 | ~ exp(-J|), as shown in Fig. 1^ with 
4; = -, = , obtained e.g. from Eq. (|7|) (for another ap- 
proach see the Appendix). At the critical point uj = uj c 
the scaling of the current in the middle of the chain de- 
pends on the parity of n. For even n it is ~ 1/n 2 , while it 
is ~ 1/n 3 for odd n. This can be theoretically seen from 



Eq.( 11 ): for odd n two quadrupoles are subtracted, lead- 



ing to a subleading-order scaling along the skew-diagonal 
of Cj.k- For uj = uj c the system is therefore also an in- 
sulator in the thermodynamic limit, however, the limit 
is reached in an algebraic way. In a way, the most in- 
teresting regime is for uj < uj c . If, when one increases 



n one "sits" at a certain resonance (thereby changing u) 
with n), ujp-m — £p + £m, and one at the same time 
also lets e — > 0, the scaling is |j(„+i)/2| ~ l/n, simply 
because only the resonance term in Eq. ([6]) contributes. 
However, the limit n — > oo and at the same time e — > 0, 
while (J depends on n, is rather artificial. More impor- 
tant is the limit when one fixes ui and e, letting only 
n — > oo. The results, shown in Fig. [3] for extremely 
large n ~ 10 9 , indicate that the average decay of the cur- 



rent is |j(„ 



+ l)/2l 



41 C, 



(n+l)/2,(n+l)/2+ll 



1A 



n, as ar- 



gued analytically (10 1. Such anomalous scaling, although 



being rather common in classical systems, has been ob- 
served only recently in a quantum setting, namely in the 
isotropic Heisenberg model |15j . 
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case for u> < lo c . In Fig.Hwe show a density plot of |C. fc | 
for u> — 8 and n — 257 as well as n = 256. One can nicely 
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FIG. 4. (Color online) Exact correlations at the non- 
equilibrium phase transition point, ui — 8. The coupling is 
e — 0.1. The C~ k is approximately given by two quadrupole 
sources at two opposite corners (high intensity in the figure) . 
For odd n (left figure) two quadrupole sources are subtracted, 
for even n (right figure) they are added. 



see very high intensity near the two opposite corners, 
reflecting the fact that the correlations can be approx- 
imately described by two quadrupole sources. For odd 
n these sources are subtracted, resulting in a suppressed 
correlations along the skew-diagonal, leading there to a 
~ 1/n 3 scaling. For even n though the quadrupoles are 
added, resulting in higher correlations scaling as ~ 1/n 2 
along the skew-diagonal. Also compare these plots with 
the cross sections along the diagonal (magnetization) and 
near-diagonal (current) shown in the Fig. ^1 
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FIG. 3. (Color online) Scaling of the current in the middle of 
the chain on n. Top left: ui > 8, top right: w = 8, bottom: 
uj — 1. Note that the scaling for the critical w, being either 
~ 1/n 2 or ~ 1/n 3 , depends on the parity of the system size. 
For uj < uj c the dependence on n is for small e not monotonic 
~ 1/y/n, however, oscillations get smaller for larger e; "Nu- 
meric" points are numerically exact solutions, "Analytic" are 
weak-coupling summation of Eq. fijj. For u) > uj c the current 
scales as ~ exp (—5*) with the evanescence length £ = 
(inset). 



VS 



IV. 



SPATIAL DEPENDENCE OF 
CORRELATIONS 

A. Critical lu c 



At the critical driving frequency u) = lo c the spatial de- 
pendence of C~ k does not show any oscillations, as is the 



B. Correlations for uj < u c 

In Fig. [5] we show the expectations of all nonzero 2- 
point observables, i.e., C~ k , at smaller frequency than u c . 
We can see that the structure of correlations can be quite 
different at different frequencies. Below the transition 
point and for small e the correlations are dominated by 
resonances, for which the condition lo = W p _ m = e p + e m 
(e p = 4 cos ^rj) is fulfilled. As explained, the position of 
resonances is different for even and odd sizes n because 
in one case only even p + m are allowed while in the other 
only odd p + m occur. This is illustrated in Fig. [6j 

If the driving frequency uj is equal to some resonant 
frequency u> p _ m , where p + m = n (mod 2), and if e is 
sufficiently small, the spatial pattern of correlations is 
given by the theoretical formula 

|C7 fc | oc I sin (pirx) sin (miry) + sin (mitx) sin (piry)\. 

This form comes due to a combination of two eigenfunc- 
tions. How small must the coupling be for this to happen 
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FIG. 5. (Color online) Correlations \C~ k \/e for n 
e = 0.001. On the left is for an out-of-resonance u — 6.163 
with an essentially "random" structure of C; on the right is 
for u) — U2-9 ~ 7.975, which is on the resonance, resulting in 
a simple form C ~ | sin (27rx) sin (9ny) + sin (9ttx) sin (2iry)\. 
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FIG. 6. Resonances are for odd and even n at different posi- 
tions. Current at the middle of the chain is for n — 256 by 
two orders of magnitude larger than for n — 257 because we 
are close to uj c = 8, where one scales as ~ l/^ 2 , while the 
other is ~ 1/n 3 . Both data are for e = 0.01. 



is a nontrivial question. It depends on how well is the res- 
onance in question resolved. For instance, the density of 
resonances is in general higher at smaller uj. This means 
that in order to resolve a resonance one will typically 
need a smaller e at smaller frequencies than for instance 
just below the critical frequency. We illustrate this phe- 
nomenon in Fig. [9] Two resonances are shown: the res- 
onance 5 — 6 at a->5_6 ~ 7.98192 is well separated from 
the others, while the 2 — 9 resonance at W2-9 ~ 7.97482 
is very close to the uje-7 ~ 7.97481. In fact, in the plot 
of |j(n+i)/2| on w these two resonances can not be re- 
solved on the scale of the plot. Because of that, at large 
e = 1 the spatial dependence of correlations is a kind 
of merger of resonances W2-9 and u>6-7- The theoretical 
small-e shape of the resonance is resolved only around 
£ = 0.001. Spatial pattern of the resonance Ws_6 is on 
the other hand well resolved already at large e = 1. 



For small ui the transition to theoretical \C~ k \ typically 
happens at smaller couplings. This is shown in Fig. [7] 
where we show the 3 — 82 resonance. This resonance is 
the closest to the frequency ui = 6.163, whose correla- 
tions have been shown in Fig. [5] We can see that one 
needs e = 0.0001 in order to reach Eq.( 12 1. The shape of 



correlations at general small frequencies and large cou- 
plings typically looks rather "random" , similar to plots 
at e = 1.0 ore = 0.01 in Fig. [7] . 
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FIG. 7. (Color online) Resonance 0J3_82 at various e. This 
resonance is the closest to w = 6.163 shown in Fig. [5] We plot 

\Cj :k \/e. 



V. OTHER RESERVOIRS 

The physical phenomenon described for our exactly 
solvable model, i.e. a phase transition with u>, is robust 
to small changes of the model. For instance, in this Sec- 
tion we show that similar transition is obtained also for 
a system coupled to the so-called two-spin baths, in each 
of which one has 16 Lindblad operators acting on two 
boundary sites (simulating finite-" temperature" ) . Phys- 
ical picture, explaining why the phase transition occurs, 
is clear. If driving is faster than the fastest time scale 
in the systems - the nearest neighbor coupling (hop- 
ping strength) - then the system is an insulator. We 



therefore conjecture that the phenomenon is robust to, 
for instance, changes of the reservoir. The only neces- 
sary ingredient is that one has reservoirs of magneti- 
zation, or, in fermionic picture, of particles. The ad- 
vantage of reservoirs studied in the present work, i.e., 
with Lindblad operators L ~ a , is that we are able 
to analytically solve for the nonequilibrium steady state. 
For other types of reservoir Lindblad operators the sys- 
tem is not quadratic in fermionic variables anymore and 
exact solution is in general not possible. Nevertheless, 
to test our conjecture we performed numerical simula- 
tions with time-dependent density matrix renormaliza- 
tion group (tDMRG) method. We used the so-called 
two-spin reservoirs, in which there are 16 Lindblad oper- 
ators acting on two spins at each chain end. Details of 
our implementation can be found in Ref. [16] . Because 
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FIG. 8. (Color online) Spin current in the middle of the chain 
for different frequencies. Symbols, squares for n = 8 and 
circles for n = 16, are results of tDMRG simulations, two 
full lines are analytic results for L ~ a . Open symbols 
are tDMRG simulation of the system studied in the present 
work (L ~ a ) while full symbols are for a two-spin bath, 
simulating finite-temperature reservoirs, e — 1, (j,q — 0.2. 



experimental realization, for instance with ions, is eas- 
ier for short chains we choose a short chain of length 
n = 8 and n = 16, in order to demonstrate that the phase 
transition can be seen already in small systems. Results 
of numerical simulations are shown in Fig. [H] First, we 
cross-checked our analytical results for Lindblad opera- 
tors L ~ a ± presented in this work with the results of 
tDMRG simulations (open symbols). By comparing two 
points at u ~ 6 and uj « 8.4 one can see a large drop in 
the current around co c = 8. Using a two-spin bath one 
can achieve a nonequilibrium steady state with nonzero 
average energy density. Assuming that the state is locally 
close to a canonical one, one can ascribe temperature to 
the nonequilibriums steady state by comparing the aver- 
age energy density with the canonical one [17] . In our 
simulations (full symbols in Fig. [8]) the energy density in 
the steady state is (cr*er* +1 + er|o-J +1 ) ~ —0.26, which 
corresponds to the canonical expectation value at tem- 
perature T ~ 7.6. The nonequilibrium steady state of 



the analytical solution has on the other hand zero aver- 
age energy density and can be described as a state at an 
infinite temperature. Observing data points for two-spin 
bath and n = 16 (full circles) at w ~ 6 and atww 8.4 we 
can see that even for a non-solvable finite-temperature 
NESS there is still a two orders of magnitude drop in the 
current as one increases the driving frequency beyond w c . 
Phase transition is therefore quite robust and is not only 
a property of our solvable model. 



VI. EXPERIMENTAL IMPLEMENTATION 

Relatively novel way of implementing various quan- 
tum models is via a rapidly developing field of simulat- 
ing physical systems with cold atoms or ions. All in- 
gredients necessary to implement our model, like the ex- 
change interaction between nearest neighbors, have al- 
ready been achieved 18 . Simulating Lindblad equation, 
for instance in order to dissipatively prepare a given pure 
state [Hfl,l2"P] . is by now also quite established. We shall 
here sketch the implementation with ions, although other 
realizations, for instance with atoms in an optical lat- 
tice [35], go along similar lines. The method would ac- 
tually be very similar to the one used in Ref. [201 [H] 
where a master equation with the Lindblad operator 
L' = \cr\(t — ofcrfcrferj) has been implemented on a 
system of 4 ions. In our model we instead need L = a\ . 
Difference from Refs. [2PJ 121] would be that to imple- 
ment L one has to apply an entangling M0lmer-S0rensen 
gate [23] just on 2 ions (ancilla and one system's site) 
instead of on 5 ions - for description see e.g. the ap- 
pendix B. in Ref. 21j. This can be seen by writing 
L = a^ = tj(Ji(1 —erf). Up to a trivial rotation around 
y- axis this is similar to V . Phase transition could be de- 
tected by simply measuring the state (magnetization) of 
ions. Another feasible way to implement our model is in 
a mesoscopic setting with electrons [21] • One would need 
a quantum wire (coupled quantum dots, molecular wire, 
etc.) coupled to electron reservoirs. By varying electro- 
chemical potential in the reservoirs, for instance by an 
external gate potential or a laser pulse, one can achieve 
a time-dependent occupation of reservoirs, which would 
correspond to our magnetization reservoirs in the spin 
language. In fact, a somewhat related model of tight- 
binding electrons in an a.c. electric field [5] has been 
studied extensively by approximate theoretical and nu- 
merical methods, see e.g. |25j and references therein. 



VII. CONCLUSION 

We have analytically solved a system of ID spinless 
fermions under harmonic a.c. driving at the lattice ends. 
With the driving frequency there is a transition from a 
system with a ballistic transport for u) = 0, to the one 
with anomalous transport for oj < 8, which at the critical 
frequency u = 8 changes to an insulator. 
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FIG. 9. (Color online) Transition to the theoretical correlation function (121 at high resonant frequencies. On the left a case 
of the 2 — 9 resonance is shown, while the right pictures show the 5 — 6 resonance. Top middle picture shows the theoretical 
resonance shape, Eq. ( |12[ ), for p — 2, m — 9. See text for more details. Color scale is adjusted between the plots for better 
detail, all is for n — 257. 
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Appendix A: Dependence on the coupling strength 



The coupling strength e between baths and the system 
is a rather trivial parameter that does not affect the main 
qualitative features of the transport. It of course influ- 
ences the size of correlations Cj t k- 



10 we show 



In Fig. 

the dependence of the current in the middle of the chain, 
ji n +i)/2 on £• The functional form is always very similar 
to the one at uj — 0, which is known exactly [SJ [6], and 
is 0'fe) = ~tt- More important is the influence of the 
value of £ on the width of resonances and whether the 
latter are isolated or not, which is in turn reflected in the 
spatial pattern of Cj k. 
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FIG. 10. Dependence of the current in the middle of the sys- 
tem on the coupling strength e. For all ui the overall functional 
dependence is similar to the one at w = (dashed curve), only 
the size of the current changes (note different prefactors for 
u/ 1.0); all is for n = 257. 



Appendix B: Exact solution of the weak coupling 
limit 



The stationary correlation matrix C 
the continuous Lyapunov equation 



is obtained from 



{A,CT} = -4£S, A = 2(J + ieR)- -1, (Bl) 

with J being an n x n matrix with the only nonzero 
elements Jk.k+i = Jk+i.k — 1 while the n x n matri- 
ces S, R have all elements equal to zero except R\^\ — 
Rn,n = Si,i = (—l) n S n ,n ~ 1- The matrix A consists 
of a term J, which corresponds to kinetic energy of the 
system, a coupling term R, that corresponds to the dissi- 
pation to the environment, and one additional term, that 
comes from the time dependent, oscillatory part of the 
correlation matrix. The source term S is determined by 
the driving, i.e. the difference in the chemical potentials 
between the left and right reservoirs. The equation (Bl) 
is straightforwardly solved by the following ansatz 



where A m = /3,V + 0W ■ 

Appendix C: Exact solution of the weak coupling 
limit for u > uj c 



In the case u> > ui c the ansatz ( B2 ) has no singularities 



and becomes analytic also in the zeroth order in e, i.e. by 
taking A m sw /3^, and for all system sizes n. Therefore, we 
can find the solution of the equation (Bl ) using a different 
approach. First we simplify the continuous Lyapunov 



equation ( B 1 1 



2{J,CT}-u;C- 



-4eS. 



(CI) 



The exact solution of the equation (CI) is sought in the 
form of a perturbative ansatz 



C = 



oo 1 



(C2) 



3=0 



3>k 



A* A 



right 



n 



an. 



(B2) 



eht 



where \&_- denotes the j-th right eigenvector of A with 



right 



ft* 



riJht 



the corresponding eigenvalue /3j, A.}&_ 

Then, by plugging the ansatz into the equation (Bl) we 

get the expression for the coefficients Aj t k 



A 



4t 



i>fe 



ft + ft 



-qj 



left* 



S* 



loft* 



(B3) 



We find the eigenvectors and the eigenvalues perturba- 
tively in e and calculate only the leading contribution to 
the correlation matrix C~ . The leading-order eigenvec- 
tors and the corresponding eigenvalues are 



The straightforward recursive solution is 



Cj — {23, Cj_i}, Co 



4eS 

U! 



(C3) 



From the recursion (C3) it is possible to calculate each 
term in the solution (C2) explicitly. After a tedious cal- 



culation we then rewrite the correlations as follows 



C U = 4e 



E 

l rn— — oo v-- 



E 



(-1)* 



G 



j—2l—i/,k—2rn,— U) 



where 



(C4) 



/(0) 



where a^ 



M 



sin ii ; i, - (3)-' 1 = — - 4cosaj, (B4) 
= 8 the sum of the eigen- 



For uj < ui c 



values in the denominator of the ansatz (B2 1 can vanish, 



hence we need to take into account the first order correc- 
tion to the eigenvalues as well 

8k 



P 



(i) 



(n + 1) 



sm" an 



(B5) 



Combining the equations (|B2|, (B3), (B4), and (B5) we 



immediately get the correlation matrix in the leading or- 
der in e, 



G jtk = jkoj-3- k T{j + k- l)T(j + k + 1) x (C5) 

16" 



4-P3 



j+k-l j+k j+fc+1 j+fc+2 
2 ' 2 ' 2 ' 2 

j + l,k + l,j + k + l 



by means of the standard Gamma function T(x) and gen- 
eralized hypergeometric function 4F3. Note that for large 
system sizes n 3> 1 and for elements of the covariance 
matrix lying near the diagonal \j — k\ <C n we can ap- 
proximate 



C 



j,k 



(G ilfc + (-l) B G i _i, fc _i). 



(C6) 



C lk = 



-32e 



E 



sm CLj p sm a& Tr 



p,m—l 

p-\-m=n(mod 2) 



(n+l) 2 (A p + A m ) 



(B6) 



From the approximate solution ( C6 1 it is possible to ex- 



tract the scaling of the elements near the diagonal in 

the limit j ^> 1. We find that Cj . +1 oc e"^'^, where 

£ = , 1 

S y/uJ — UJ C 
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